

#include "glodef.h"
#include "errlog.h"
#include "grid3d.h"
#include "g3dutils.h"
#include "vecs.h"
#include "fe.h"
#include "laplace.h"
#include "f3_parser.h"
#include "alegrid.h"

#include <stdlib.h>
#include <string.h>

static inline void _ale_error_alloc(const char *func)
{
    error_msg("alegrid: memory allocation error in `%s'\n", func);
}

/****************************************************************************
 * Create an ale grid and return it
 */
alegrid_t *init_alegrid(grid3d_t *grid)
{
    alegrid_t *ag;

    PrintLog("Initializing alegrid.");
    ag=calloc(1, sizeof(alegrid_t));
    if(ag==NULL)
        _ale_error_alloc(__func__);

    if(grid==NULL)
        error_msg("Must supply a valid grid when calling `%s'.\n", __func__);
    else
        ag->grid=grid;

    /* Allocate memory for members of ag */
    ag->ref_nodes=calloc(grid->NN, sizeof(node_t));
    if(ag->ref_nodes==NULL)
        _ale_error_alloc(__func__);

    ag->old_nodes=calloc(grid->NN, sizeof(node_t));
    if(ag->old_nodes==NULL)
        _ale_error_alloc(__func__);

    ag->ndisp=alloc_vec(grid->NN);

    ag->can_move=malloc(sizeof(int)*grid->NN);
    if(ag->can_move==NULL)
        _ale_error_alloc(__func__);

    /* Save current nodes positions in ref_nodes */
    memcpy(ag->ref_nodes, grid->nodes, grid->NN*sizeof(node_t));
    /* all nodes on the boundary have BC_DIR when solving Laplace equation
     * for the moving grid.
     */
    {
        int i;
        for(i=0; i<grid->NN; i++)
            if(ag->ref_nodes[i].on_bdry)
            {
                ag->ref_nodes[i].bcond=BC_DIR;
                ag->can_move[i]= -1;
            }
            else
                ag->can_move[i]=1;
    }

    PrintLog(" done.\n");
    /* initialize the Laplace solver */
    {
        unsigned eltype=fe_p1;
        node_t *nodes=grid->nodes;
        grid->nodes=ag->ref_nodes; /* same coords, but bcond may be different */
        UpdateMeasures(grid);
        InitFE(1, &eltype, grid);
        ag->laplace_data=InitLaplace(grid, eltype);
        /* restore grid->nodes */
        grid->nodes=nodes;
        UpdateMeasures(grid);
    }

    return ag;
}

/****************************************************************************
 * Routines that change the ref_nodes array
 */

void set_ref_nodes(alegrid_t *ag, node_t *new_ref_nodes)
{
    grid3d_t *grid=ag->grid;
    node_t *arn=ag->ref_nodes;
    for(int i=0; i<grid->NN; i++)
    {
        memcpy(arn->coords, new_ref_nodes->coords, 3*sizeof(double));
        ++arn;
        ++new_ref_nodes;
    }
    /* save grid->nodes */
    node_t *save_nodes=grid->nodes;
    grid->nodes=ag->ref_nodes;
    UpdateMeasures(grid);
    /* reassemble matrix */
    ReassembleLaplace(ag->laplace_data);
    /* restore grid->nodes */
    grid->nodes=save_nodes;
    UpdateMeasures(grid);
}


/****************************************************************************
 * Assign boundary conditions for the moving parts of the boundary
 */

#define ALEBC_ID        0x68DE1B1A

typedef struct alebc_s
{
    unsigned struct_id;
    struct alebc_s *next;
    unsigned on_bdry;
    char * ptpos_s;
    field3d *ptpos_r;
} alebc_t;

/****************************************************************************
 * Release an ale grid
 */
void done_alegrid(alegrid_t *alegrid)
{
    if(alegrid==NULL)
        return;
    free(alegrid->can_move);
    free(alegrid->ref_nodes);
    free(alegrid->old_nodes);
    free_vec(alegrid->ndisp);
    {
        alebc_t *abc=alegrid->alebc_data;
        while(abc)
        {
            alebc_t *a=abc->next;
            free(abc);
            abc=a;
        }
    }
    DoneLaplace(alegrid->laplace_data);
    free(alegrid);
    PrintLog("alegrid done.\n");
}

vector get_ndisp(alegrid_t *ag)
{
    vector nd;
    node_t *cn, *on;
    int i;
    if(ag==NULL)
        return NULL;
    assert(ag->grid!=NULL);
    assert(ag->old_nodes!=NULL);
    assert(ag->ndisp!=NULL&&get_N_vec(ag->ndisp)==ag->grid->NN);

    cn=ag->grid->nodes;
    on=ag->old_nodes;
    nd=ag->ndisp;

    /*
        if(mpi_my_id == 0)
        {
            for(i=0;i<get_N_vec(nd);i++)
            {
                double vals[3];
                for(int j=0;j<3;j++)
                    vals[j] = cn[i].coords[j] - on[i].coords[j];
                set_value_vec(nd, i, vals);
            }
        }
     */

    use_vector(nd)
    for(i=0; i<get_N_vec(nd); i++)
    {
        ndX[i]=cn[i].coords[0]-on[i].coords[0];
        ndY[i]=cn[i].coords[1]-on[i].coords[1];
        ndZ[i]=cn[i].coords[2]-on[i].coords[2];
    }
    done_vector(nd)

    return nd;
}

static alebc_t *_get_alebc(void **r, unsigned on_bdry)
{
    alebc_t *root= *r;

    if(root==NULL)
    {
        root=calloc(1, sizeof(alebc_t));
        if(root==NULL)
            _ale_error_alloc(__func__);
        root->struct_id=ALEBC_ID;
        root->on_bdry=on_bdry;
        (*r)=root;
        return root;
    }
    else
    {
        alebc_t *abc=root;

        if(root->struct_id!=ALEBC_ID)
            error_msg("alegrid: `r` is not a pointer to an "
                      "`alebc_t` structure.\n");

        while(abc->next)
        {
            if(abc->on_bdry==on_bdry)
                return abc;
            abc=abc->next;
        }
        if(abc->on_bdry==on_bdry)
            return abc;
        abc->next=calloc(1, sizeof(alebc_t));
        if(abc->next==NULL)
            _ale_error_alloc(__func__);
        abc=abc->next;
        abc->struct_id=ALEBC_ID;
        abc->on_bdry=on_bdry;
        return abc;
    }
}

/*****************************************************************************
 *
 */
static void compute_alebc(const alegrid_t *ag, unsigned on_bdry, double T,
                          const double *Pt, double *V)
{
    alebc_t *abc;

    if(ag==NULL)
        error_msg("alegrid is NULL in `%s'\n", __func__);
    memcpy(V, Pt, 3*sizeof(double));
    abc=ag->alebc_data;
    while(abc!=NULL)
    {
        if(abc->struct_id!=ALEBC_ID)
            error_msg("Data corruption 1 detected in `%s'\n", __func__);
        if(abc->on_bdry&on_bdry)
        {
            if(abc->ptpos_r)
            {
                abc->ptpos_r(T, Pt, V);
                return;
            }
            else if(abc->ptpos_s)
            {
                if(!parse_vect(T, Pt, V, abc->ptpos_s))
                    error_msg("Failed to parse alebc string in `%s'\n", __func__);
                return;
            }
            else
                error_msg("Data corruption 2 detected in `%s'\n", __func__);
        }
        abc=abc->next;
    }
}

/****************************************************************************
 * Move grid nodes by solving Laplace equation
 */

/* main driver routine */
static void __move_alegrid_laplace_bc_str
(alegrid_t *ag, double time,
 void (*bc_func)(alegrid_t *ag, vector rhs, vector newpos, double time))
{
    int NN, i;
    vector rhs, newpos, tmp;
    node_t *nodes;
    grid3d_t *grid;

    if(ag==NULL)
        error_msg("alegrid is NULL in `%s'\n", __func__);
    if(ag->grid==NULL)
        error_msg("grid is NULL in `%s'\n", __func__);
    grid=ag->grid;
    NN=grid->NN;

    rhs=alloc_vec(NN);
    tmp=alloc_vec(NN);
    newpos=alloc_vec(NN);

    /* copy current nodes' positions into old_nodes */
    /* NOTE: must copy since this copies boundary codes as well. If we just
     swap pointers, this info may become inacurate */
    memcpy(ag->old_nodes, grid->nodes, NN*sizeof(node_t));
    /* save a copy of the grid->nodes pointer */
    nodes=grid->nodes;
    /* place reference nodes in grid */
    grid->nodes=ag->ref_nodes;
    UpdateMeasures(grid);

    bc_func(ag, rhs, newpos, time);

    SolveLaplace_vector(ag->laplace_data, newpos, rhs, tmp, TRUE, FALSE);
    /* TODO : !!!!!! _CHANGE_reassemble_TO_TRUE_ABOVE_ !!!!!!!  */

    use_vector(newpos)
    for(i=0; i<NN; i++)
    {
        nodes[i].coords[0]=newposX[i];
        nodes[i].coords[1]=newposY[i];
        nodes[i].coords[2]=newposZ[i];
    }
    done_vector(newpos)

    grid->nodes=nodes;
    UpdateMeasures(grid);

    free_vec(newpos);
    free_vec(tmp);
    free_vec(rhs);
    return;
}

/****************************************************************************/

/*   positions of boundary nodes are prescribed by routines */
static void __assign_bc_str(alegrid_t *ag, vector rhs, vector newpos,
                            double time)
{
    grid3d_t *grid=ag->grid;
    int i, NN=grid->NN;

    clean_vec(rhs);

    use_vector(rhs)
    use_vector(newpos)
    for(i=0; i<NN; i++)
    {
        node_t *n=grid->nodes+i;
        /* assign initial condition */
        newposX[i]=n->coords[0];
        newposY[i]=n->coords[1];
        newposZ[i]=n->coords[2];
        /* assign boundary conditions */
        if(n->bcond==BC_DIR)
        {
            double val[3];
            compute_alebc(ag, n->on_bdry, time, n->coords, val);
            newposX[i]=rhsX[i]=val[0];
            newposY[i]=rhsY[i]=val[1];
            newposZ[i]=rhsZ[i]=val[2];
        }
    }
    done_vector(newpos)
    done_vector(rhs)
}

void move_alegrid_laplace(alegrid_t *ag, double time)
{
    __move_alegrid_laplace_bc_str(ag, time, __assign_bc_str);
}

/****************************************************************************/
/*   positions of boundary nodes are given in an input vector */
static vector __bnpos;

static void __assign_bc_bnpos(alegrid_t *ag, vector rhs, vector newpos,
                              double time)
{
    grid3d_t *grid=ag->grid;
    int i, NN=grid->NN;

    assert(__bnpos!=NULL&&get_N_vec(__bnpos)==get_N_vec(rhs));

    clean_vec(rhs);

    use_vector(rhs)
    use_vector(newpos)
    use_vector(__bnpos)
    for(i=0; i<NN; i++)
    {
        node_t *n=grid->nodes+i;
        /* assign initial condition */
        newposX[i]=n->coords[0];
        newposY[i]=n->coords[1];
        newposZ[i]=n->coords[2];
        /* assign boundary conditions */
        /* assign boundary conditions */
        if(n->bcond==BC_DIR)
        {
            newposX[i]=rhsX[i]=__bnposX[i];
            newposY[i]=rhsY[i]=__bnposY[i];
            newposZ[i]=rhsZ[i]=__bnposZ[i];
        }
    }
    done_vector(__bnpos)
    done_vector(newpos)
    done_vector(rhs)
}

void move_alegrid_laplace_bnpos(alegrid_t *ag, double time, vector bnpos)
{
    __bnpos=bnpos;
    __move_alegrid_laplace_bc_str(ag, time, __assign_bc_bnpos);
    __bnpos=NULL;
}

void SetALEBCString(alegrid_t *ag, unsigned on_bdry, const char *ptpos_s)
{
    alebc_t *abc;
    if(ag==NULL)
        error_msg("alegrid is NULL in `%s'\n", __func__);
    abc=_get_alebc(&ag->alebc_data, on_bdry);
    abc->ptpos_r=NULL;
    if(abc->ptpos_s)
        free(abc->ptpos_s);
    if(ptpos_s==NULL||ptpos_s[0]=='\0')
        abc->ptpos_s=NULL;
    else
    {
        int len=strlen(ptpos_s)+1;
        char *s=(char *) calloc(len, sizeof(char));
        if(s==NULL)
            _ale_error_alloc(__func__);
        abc->ptpos_s=strncpy(s, ptpos_s, len);
    }
}

void SetALEBCRoutine(alegrid_t *ag, unsigned on_bdry, field3d *ptpos_r)
{
    alebc_t *abc;
    if(ag==NULL)
        error_msg("alegrid is NULL in `%s'\n", __func__);
    abc=_get_alebc(&ag->alebc_data, on_bdry);
    abc->ptpos_r=ptpos_r;
    if(abc->ptpos_s)
    {
        free(abc->ptpos_s);
        abc->ptpos_s=NULL;
    }
}







